Objetivo

Fonte dos dados:

Carregar pacotes

library(tidyverse)
library(knitr)
library(lubridate)
library(janitor)
library(broom)
library(gghighlight)
library(skimr)
library(DT)

Carregar a base de dados

# ler o arquivo
mananciais <-
  read_delim(
    "https://github.com/beatrizmilz/mananciais/raw/master/inst/extdata/mananciais.csv",
    delim = ";",
    escape_double = FALSE,
    col_types = cols(data = col_date(format = "%Y-%m-%d")),
    locale = locale(decimal_mark = ",", grouping_mark = "."),
    trim_ws = TRUE
  )

Conhecendo a base de dados

Quais colunas estão disponíveis?

glimpse(mananciais)
## Rows: 50,997
## Columns: 8
## $ data                <date> 2022-07-12, 2022-07-12, 2022-07-12, 2022-07-12, 2…
## $ sistema             <chr> "Cantareira", "Alto Tietê", "Guarapiranga", "Cotia…
## $ volume_porcentagem  <dbl> 38.4, 58.1, 72.1, 79.0, 95.5, 43.0, 85.3, 38.5, 58…
## $ volume_variacao     <dbl> -0.1, -0.2, -0.4, -0.2, -0.1, -0.2, -0.3, -0.1, -0…
## $ volume_operacional  <dbl> 376.97898, 325.41513, 123.36649, 13.03291, 107.103…
## $ pluviometria_dia    <dbl> 0.0, 0.0, 0.0, 0.2, 0.2, 0.2, 0.2, 0.1, 0.2, 0.0, …
## $ pluviometria_mensal <dbl> 0.3, 1.9, 0.2, 1.2, 0.6, 2.8, 1.2, 0.3, 1.9, 0.2, …
## $ pluviometria_hist   <dbl> 46.5, 46.9, 41.5, 51.3, 54.2, 91.3, 77.7, 46.5, 46…

Sumário da base de dados

skim(mananciais)
Data summary
Name mananciais
Number of rows 50997
Number of columns 8
_______________________
Column type frequency:
character 1
Date 1
numeric 6
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
sistema 0 1 5 12 0 7 0

Variable type: Date

skim_variable n_missing complete_rate min max median n_unique
data 0 1 2000-01-01 2022-07-12 2011-08-20 8229

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
volume_porcentagem 0 1 66.57 25.72 -24.30 47.70 69.20 88.7 110.50 ▁▂▇▇▇
volume_variacao 0 1 0.00 0.59 -6.40 -0.20 -0.10 0.1 16.60 ▁▇▁▁▁
volume_operacional 0 1 151.27 185.47 -238.84 13.96 97.02 200.8 987.05 ▁▇▂▁▁
pluviometria_dia 0 1 4.05 9.86 0.00 0.00 0.20 2.8 222.50 ▇▁▁▁▁
pluviometria_mensal 0 1 64.14 71.03 0.00 9.80 41.20 94.4 557.60 ▇▂▁▁▁
pluviometria_hist 245 1 132.26 68.40 23.50 75.00 123.40 189.1 298.90 ▇▇▅▅▂

Alterando a base de dados

Criar a coluna ano, mes_ano, e ordenar a coluna de sistema:

mananciais <- mananciais |>
  mutate(
    ano = year(data),
    mes_ano = floor_date(data, "month"),
    sistema = fct_relevel(
      sistema,
      levels = c(
        "Cantareira",
        "Alto Tietê",
        "Guarapiranga",
        "Rio Grande"  ,
        "São Lourenço" ,
        "Cotia" ,
        "Rio Claro"
      )
    )
  )

Explorando os dados:

Quais são os sistemas existentes na base? Quantas observações temos para cada Sistema?

mananciais |> 
  count(sistema) |> 
  kable()
sistema n
Cantareira 8229
Alto Tietê 8229
Guarapiranga 8229
Rio Grande 8229
São Lourenço 1623
Cotia 8229
Rio Claro 8229 
mananciais |> 
  group_by(sistema) |> 
  filter(data == min(data)) |> 
  select(sistema, data) |> 
  kable()
sistema data
São Lourenço 2018-02-01
Cantareira 2000-01-01
Alto Tietê 2000-01-01
Guarapiranga 2000-01-01
Cotia 2000-01-01
Rio Grande 2000-01-01
Rio Claro 2000-01-01

Quais são os sistemas com o maior volume operacional armazenado?

mananciais |> 
  filter(data == max(data)) |> 
  arrange(desc(volume_operacional)) |> 
  select(sistema, volume_operacional) |> 
  kable()
sistema volume_operacional
Cantareira 376.97898
Alto Tietê 325.41513
Guarapiranga 123.36649
Rio Grande 107.10315
São Lourenço 75.72723
Cotia 13.03291
Rio Claro 5.87214

Gráficos!

Gráficos de linhas

mananciais |> 
  filter(sistema == "Cantareira") |> 
  ggplot() +
  geom_line(aes(x = data, y = volume_porcentagem)) +
  theme_minimal() +
  labs(x = "Ano", y = "Volume (%)")

mananciais |>
  filter(sistema == "Cantareira") |>
  ggplot() +
  geom_line(aes(x = data, y = volume_porcentagem)) +
  theme_minimal() +
  labs(x = "Ano", y = "Volume (%)") +
  gghighlight(ano %in% c(2014:2015))

  # gghighlight(mes_ano >= "2013-04-01", mes_ano < "2016-01-01")

Exercício 1

Adapte o código abaixo para criar um gráfico para o sistema Alto Tietê:

mananciais |> 
  filter(sistema == "Cantareira") |> 
  ggplot() +
  geom_line(aes(x = data, y = volume_porcentagem)) +
  theme_minimal() +
  labs(x = "Ano", y = "Volume (%)")

Criando uma tabela de resumo!

mananciais |> 
  summarise(min = min(volume_porcentagem),
            max = max(volume_porcentagem),
            media = round(mean(volume_porcentagem), 1),
            variancia = round(var(volume_porcentagem), 1), 
            desvio_padrao = round(sd(volume_porcentagem), 1)) |> 
  kable()
min max media variancia desvio_padrao
-24.3 110.5 66.6 661.8 25.7 
mananciais |> 
  group_by(ano, sistema) |> 
  summarise(min = min(volume_porcentagem),
            max = max(volume_porcentagem),
            media = round(mean(volume_porcentagem), 1),
            variancia = round(var(volume_porcentagem), 1), 
            desvio_padrao = round(sd(volume_porcentagem), 1)) |> 
  datatable()

Exercício 2

Adapte o código a seguir para criar uma tabela de resumo para a variável de Volume Operacional:

mananciais |> 
  group_by(ano, sistema) |> 
  summarise(min = min(volume_porcentagem),
            max = max(volume_porcentagem),
            media = round(mean(volume_porcentagem), 1),
            variancia = round(var(volume_porcentagem), 1), 
            desvio_padrao = round(sd(volume_porcentagem), 1)) |> 
  datatable()

Regressão linear

modelo_linear <- lm(volume_variacao ~ pluviometria_dia, data = mananciais)

modelo_linear
## 
## Call:
## lm(formula = volume_variacao ~ pluviometria_dia, data = mananciais)
## 
## Coefficients:
##      (Intercept)  pluviometria_dia  
##         -0.13892           0.03504
tidy(modelo_linear) |> kable()
term estimate std.error statistic p.value
(Intercept) -0.1389243 0.0022709 -61.17589 0
pluviometria_dia 0.0350441 0.0002131 164.41319 0 
glance(modelo_linear)  |> kable()
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual nobs
0.3464416 0.3464288 0.4743759 27031.7 0 1 -34329.33 68664.66 68691.18 11475.53 50995 50997 
mananciais |>
  group_by(sistema) |>
  summarise(
    modelo = lm(volume_variacao ~ pluviometria_dia, data = cur_data()) |>
      tidy() |>
      select(term, estimate) |>
      pivot_wider(names_from = term, values_from = estimate),
    r2 = lm(volume_variacao ~ pluviometria_dia, data = cur_data()) |>
      glance() |>
      pull(adj.r.squared) |>
      round(2)
  ) |>
  unnest(cols = c(modelo)) |>
  clean_names() |>
  mutate(
    modelo = str_glue(
      "volume_variacao = {round(intercept,3)} + {round(pluviometria_dia,3)} * pluviometria_dia"
    )
  ) |>
  select(sistema, modelo, r2) |>
  arrange(desc(r2)) |>
  kable()
sistema modelo r2
Rio Claro volume_variacao = -0.279 + 0.048 * pluviometria_dia 0.51
Cotia volume_variacao = -0.125 + 0.04 * pluviometria_dia 0.41
Rio Grande volume_variacao = -0.142 + 0.035 * pluviometria_dia 0.38
Guarapiranga volume_variacao = -0.105 + 0.031 * pluviometria_dia 0.29
Alto Tietê volume_variacao = -0.074 + 0.02 * pluviometria_dia 0.28
Cantareira volume_variacao = -0.07 + 0.018 * pluviometria_dia 0.21
São Lourenço volume_variacao = -0.103 + 0.035 * pluviometria_dia 0.12 
mananciais_r2 <- mananciais |>
  group_by(sistema) |>
  summarise(
    r2 = lm(volume_variacao ~ pluviometria_dia, data = cur_data()) |>
      glance() |>
      pull(adj.r.squared) |>
      round(2),
    max_pluviometria_dia = max(pluviometria_dia),
    max_volume_variacao = max(volume_variacao),
  )


mananciais |>
  ggplot(aes(y = volume_variacao, x = pluviometria_dia)) +
  geom_point(aes()) +
  geom_smooth(method = "lm") +
  facet_wrap(vars(sistema)) +
  geom_text(aes(
    x = max(max_pluviometria_dia)*0.80,
    y = max(max_volume_variacao)*0.9,
    label = paste0("R2 = ", r2)
  ), data = mananciais_r2) +
  theme_minimal() +
  labs(x = "Pluviometria do dia", y = "Variação do volume diário")

Exercício 3

manancias_filtrado <- mananciais |>
  filter(ano >= 2000)

mananciais_r2_filtrado <- manancias_filtrado |>
  group_by(sistema) |>
  summarise(
    r2 = lm(volume_variacao ~ pluviometria_dia, data = cur_data()) |>
      glance() |>
      pull(adj.r.squared) |>
      round(2),
    max_pluviometria_dia = max(pluviometria_dia),
    max_volume_variacao = max(volume_variacao),
  )


manancias_filtrado |> 
  ggplot(aes(y = volume_variacao, x = pluviometria_dia)) +
  geom_point(aes()) +
  geom_smooth(method = "lm") +
  facet_wrap(vars(sistema)) +
  geom_text(aes(
    x = max(max_pluviometria_dia)*0.80,
    y = max(max_volume_variacao)*0.9,
    label = paste0("R2 = ", r2)
  ), data = mananciais_r2_filtrado) +
  theme_minimal() +
  labs(x = "Pluviometria do dia", y = "Variação do volume diário")

Para saber mais